Non-monotonic temperature evolution of dynamic correlations in glass-forming liquids 
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The viscosity of glass-forming liquids increases by many orders of magnitude if their temperature 
is lowered by a mere factor of 2-3^'^ . Recent studies suggest that this widespread phenomenon 
is accompanied by spatially heterogeneous dynamics^, and a growing dynamic correlation length 
quantifying the extent of correlated particle motion^"^. Here we use a novel numerical method 
to detect and quantify spatial correlations which reveal a surprising non-monotonic temperature 
evolution of spatial dynamical correlations, accompanied by a second length scale that grows mono- 
tonically and has a very different nature. Our results directly unveil a dramatic qualitative change in 
atomic motions near the mode-coupling crossover temperature^which involves no fitting or indirect 
theoretical interpretation. Our results impose severe new constraints on the theoretical description 
of the glass transition, and open several research perspectives, in particular for experiments, to 
confirm and quantify our observations in real materials. 
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More than forty years ago Adam and Gibbs^ put for- 
ward the seminal idea that the relaxation dynamics of 
highly viscous liquids occurs through 'cooperatively re- 
laxing regions'. Their theory suggested that particle mo- 
tion occurs in a collective manner in localized domains^ 
whose typical size is related to the entropy of the systems 
and increases with decreasing temperature. This result 
implied that the relaxation mechanism is controlled by a 
unique length scale of thermodynamic origin with a di- 
rect signature in the dynamics. This approach regained 
momentum in the 1990's when novel experimental tech- 
niques and large scale computer simulations established 
the presence of dynamical heterogeneities, i.e. localized 
regions where dynamics is significantly different from the 
average^ii, although these observations can also be inter- 
preted as a purely dynamical phenomenoniS. 

A qualitatively similar, but much more detailed, the- 
oretical description is obtained within the framework 
of the random first order transition (RFOT) theoryii. 
Within this approach, there exists an ideal glass transi- 
tion that underlies glass formation, with an associated 
diverging correlation length scale of entropic origin since 
it is related to the existence of a large number of long- 
lived metastable states. At very low temperatures, the 
glass-former is described as a 'mosaic' of correlated do- 
mains that rearrange in a thermally activated, collective 
manner, such that again static and dynamic correlations 
coincide and grow with the viscosityi^. However, mod- 
els with an RFOT (mainly mean-field like models) also 
display a 'spinodal' singularity at a higher temperature, 
Tci at which metastability is lost and therefore the mosaic 
picture is no longer usefuliiii^. Thus above Tc, a different 
approach must be used and it is found that the dynamics 
within mean- field models ^^ has profound similarities with 
the one predicted from mode-coupling theory^, where Tc 
corresponds to a dynamic critical point associated with 
a diverging dynamical correlation length!^. Thus, even 
at the theoretical level, the physics around Tc remains 
'mysterious'—: How can a dynamical correlation length 



diverge at two distinct temperatures? The mystery thick- 
ens in finite dimensions for which the mode-coupling sin- 
gularity is cut off, and its existence can be inferred only 
from fitting relaxation datai^', a procedure that is prone 
to criticisms, and therefore the physical relevance of the 
'avoided' singularity at Tc has remained a debated issue. 

Progress is also slow because experiments on molecular 
systems do not have enough resolution to follow atomic 
motions over long times'^, and numerical simulations of- 
ten cannot access low enough temperatures to make def- 
inite statements on the size and nature of dynamic het- 
erogeneities (in fact there are so far no numerical studies 
on the dynamical heterogeneities below Tc). At present, 
the most direct measurements^ seem to indicate that the 
dynamical correlation length increases from 1-2 particle 
diameters at moderately supercooled temperatures to 5- 
10 diameters close to the glass transition temperature, 
but the interpretation of the experimental data is often 
rather difficult^. In this article we show that the use of a 
new methodology to characterize bulk dynamic correla- 
tions reveals details on their temperature evolution not 
observed in previous work. 

In parallel to the quest for a dynamical length scale, 
evidence has also been found for an increasing of static 
correlations. However, this information is not captured 
by standard two-point correlation functions. Recent 
work has for instance suggested the growth of locally fa- 
vored geometric structures in some model systemsiii^, 
but these methods are not easily generalized to different 
glass-formers. One possibility to tackle this problem are 
the point-to-set correlations which are an elegant, gen- 
eral method to capture the multi-point static correlations 
which might characterize the non-trivial structure of vis- 
cous liquids^^. The conceptual idea is to 'pin' the posi- 
tion of a number of particles (the 'set') in an equilibrated 
configuration of the fluid, and to measure how the posi- 
tion of the remaining particles is affected. It has recently 
been argued that in the geometry in which particles out- 
side a spherical cavity are pinned, this point-to-set corre- 



lation should detect the typical domain size of the RFOT 
mosaic state^-. Numerical simulations confirmed qualita- 
tively the growth of point-to-set correlations in this par- 
ticular geometry ^"'^. However, the connection to dynamic 
correlations and the precise temperature dependence of 
the various length scales were not studied, and these re- 
sults did not resolve the 'mystery'i£ of the Tc crossover. 

Inspired by previous work on confined fiuids^^ (see also 
Ref.^'^), we have generalized the idea of a point-to-set 
correlation to a novel geometry. We pin particles in a 
semi-infinite space and detect the resulting effect on the 
other half space. The principal advantage is that we can 
measure simultaneously the static and dynamic profiles 
induced by the frozen wall. Additionally we are able 
to perform simulations from the high temperature liquid 
down to and below Tc with a realistic molecular dynam- 
ics, thus allowing us to resolve at once multi-point static 
and dynamic correlations in a very broad temperature 
regime encompassing the (hypothetical) mode-coupling 
crossover. 

We study a binary mixture of quasi-hard spheres^*, 
as described in the Methods. The fiuid is equilibrated 
by means of standard molecular dynamics using periodic 
boundary conditions in all three directions. To simulate 
particles pinned within a semi-infinite space, z < 0, it is 
enough to freeze at an arbitrary time i = the position 
of all particles within a slice of thickness dwaii = lAa 
which is perpendicular to the z— axis: They form our 
'set'. Because we use periodic boundary conditions, we 
work with a very large system size in the z-direction, 
ensuring that bulk behavior is recovered at the center 
of the simulation box, i.e. the replicated walls do not 
interfere with each other. 

To properly measure a point-to-set correlation it is cru- 
cial that the frozen walls have the same structure as the 
equilibrated liquid at temperature T, such that the av- 
erage static properties of the confined liquid are unper- 
turbed^^'^^ (see also Supplementary Material). We have 
measured how far (in z) the wall influences the static lo- 
cal density fleld and its dynamics, giving us independent 
access to static and dynamic correlation length scales. 
Note that within RFOT theory, the spatial extent of the 
static profile near a wall is not directly controlled by the 
mosaic length scalei^, and a comparison with results ob- 
tained with a spherical cavity2i is not straightforward. 

A convenient observable to characterize the influence of 
the wall is the overlap profile qdt, z), defined as foUowa^i. 
We discretize space into small cubic boxes of linear size 
5 « 0.55ct, and define ni{t) = 1 if box i is occupied 
by at least one particle at time t, and ni{t) = if not. 
The overlap profile in the z direction with respect to the 
template configuration at time i = is 



qc{t,z) 
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where the sums run over all boxes at distance z from the 
wall, (• • • ) is the usual thermal average, and [■ • -Jwaii is 



FIG. 1: (a) Time dependence of the overlap qc{t,z), Eq. ([T]), 
for T = 8.0 and various values of z, increasing from z — 0.2736 
(near the wall, rightmost curve) by increments of 0.2. The 
leftmost curve presents the bulk behavior, (b) 2— dependence 
of the nontrivial part of the static overlap qoo{z) for differ- 
ent temperatures (symbols). The solid lines are fits with the 
exponential given in Eq. ((3]). 



an additional average over independent wall realizations. 
Thus qc{t,z) quantifies the similarity of particle config- 
urations separated by a time t at distance z, and, by 
construction, qc{t = 0,z) = 1, for all z. We have also 
studied qs{t, z), the single particle version of Eq. ([1]), ob- 
tained by requesting that the box i is occupied at times 
and t by the same particle. We find no relevant dif- 
ference between these two correlation functions as far as 
time dependence is concerned. 

In Fig. [1^ we show the time dependence of qc(t, z) 
at T — 8.0 for different values of z. For large z (left- 
most curves) the correlators become independent of z 
and present the bulk behavior. The presence of a shoul- 
der at intermediate times reflects the usual cage motion 
of particles observed in glassy systems-^. In the long-time 
limit, t -^ 00, the correlator decays to grand = 0.110595, 
the probability that a box is occupied, a quantity we mea- 
sured with high precision from bulk simulations. With 
decreasing z, the height of the plateau at intermediate 
times increases and the timescale to reach it decreases. 
The final decay is much slower near the wall than in the 



bulk, showing that the a-relaxation is strongly affected 
by the frozen wall. Furthermore the long-time limit of the 
overlap increases from its trivial value, qdt — s- oo,z) > 
9rand, showing that sufhciently close to the wall, the den- 
sity cannot freely fluctuate. (See Supplementary Mate- 
rial for an illustration.) Thus the set of frozen particles 
with z < influences the position of the liquid particles 
at 2; > over a non-trivial static length scale. Finally, 
Fig. [1^ shows that there exists a range of z values for 
which the long-time limit of the overlap is the trivial 
bulk value, while the relaxation timescale is slower than 
the bulk. This directly shows, with no further analysis, 
that dynamic correlations have a larger range than static 
ones, as we confirm below. 

To quantify these observations, we fit the final decay 
of qc{t, z) to the stretched exponential form 



qdt, z) - grand = A exp[- (t/r)'^] + qo 



(2) 



where A, r, /3, and q^o are fitted for each z. The profile of 
the static overlap, (700(2) measures how far from the wall 
density fiuctuations are correlated, while r(z) measures 
how far dynamics is affected. We find that Eq. ([2]) also 
describes well the single particle overlap, qs {t, z) , with 
the obvious difference that qs{t — t- 00, z) = 0, because 
eventually all particles leave the box that they occupy at 
t = 0. Thus, we obtain a second, 'self dynamic profile 
from the study of qsit, z). 

In Fig.[T]3, we display the temperature evolution of the 
static overlap profiles qao{z)- The semi- log plot suggests 
to describe these data using an exponential decay^^ 



9eo(2) = Sexp(-z/f*-'), 



(3) 



which allows us to define a static point-to-set correlation 
length scale £,^^^^{T). From these data it is clear that £,^^^^ 
grows when temperature decreases, a result in good qual- 
itative agreement with previous work^i using a very dif- 
ferent geometry. Notice that B{T) also changes rapidly 



with T, which suggests to define ^' 



■stat — int 



B{T)^^ 



a convenient estimate of the integrated profile, ^*'*^' '"' ss 
J^ qoc,iz)dz. 

We now analyze the dynamic profiles. To take into ac- 
count the fact that the amplitude and stretching of the 
time dependent correlations evolve with z, see Fig.[T^, we 
have calculated the area under the correlators qdt, z) and 
qs{t, z), taking into account only the secondary, slowly re- 
laxing part. We denote the resulting times by Tc{z) and 
Ts{z), respectively Previous studies^ have suggested 
that for large z the z— dependence of Ts{z) can be de- 
scribed well by an exponential functional form. 



log(r,) = log(T^ 



bulk\ 



S,exp(-z/e") 



(4) 



where Bs{z) and, more importantly, the dynamic length 
scale, ^s^'^iT), are adjusted for each z. The bulk re- 
laxation time T^^^^{T) is measured independently with a 
very good precision, see Fig. [2^. Using a power law fit in- 
spired by mode-coupling theory^ t^''^^{T) - [T-T^)-^, 
we obtain Tc ~ 5.2, but deviations from the algebraic fit 
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FIG. 2: (a) Temperature evolution of bulk relaxation time 
for large (top) and small (bottom) particles. Solid lines rep- 
resent a power-law fit with Tc = 5.2 (dashed line), the arrow 
indicates T = 6.0 where deviations from the fit appear, (b) 
z— dependence of the relaxation time Ts{z, T) at various tem- 
peratures. Note the non-monotonic T— dependence at inter- 
mediate and large z. 



appear above Tc near T — 6.0, see Fig. [5^ (see Supple- 
mentary Material for a discussion of this fit.) 

In Fig. ^\d we show that the data at large z can indeed 
be fitted well by the Ansatz given by Eq. ^. We have 
found a similar behavior for Tc{z), providing us with a 
second dynamic correlation length scale, S,^^'^- For high 
temperatures Eq. Qj gives a good description of the data 
over the entire range of distances z. For intermediate and 
low temperatures we see the development of a curvature 
in a semi- log plot, indicating at small z deviations from 
the simple exponential dependence. This might suggest, 
although very indirectly, the appearance of more than 
one relaxation process for the relaxation dynamics. 

A remarkable behavior occurs at intermediate and 
large distances, which has, to our knowledge, remained 
undetected. The dynamic profiles exhibit a striking non- 
monotonic evolution with temperature. A close inspec- 
tion of Fig. [2}d shows that the dynamic profiles extend to 
increasingly larger distances when temperature decreases 
from T = 30.0 down to T = 6.0, but become shorter- 
ranged when T is decreased further, down to T = 5.5, 
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FIG. 3: Temperature dependence of static and dynamic 
length scales identified in this work. The arrow at T = 6.0 
and the dashed line at Tc. — 5.2 are as in Fig. [2] Dynamic 
length scales display a non-monotonic behavior with a maxi- 
mum at T = 6.0, while static length scales increase modestly 
above T^. 



T = 5.25 and then T = 5.0. The maximum occurs near 
T — 6.0, which is also the temperature at which devia- 
tions from mode-coupling fits appear, see Fig. [2^. 

We have carefully checked that this behavior is not a 
result of our numerical analysis. A direct visual inspec- 
tion of the time correlations functions Qs (z, t) reveals that 
the spread of the curves in the slow decay has a maxi- 
mum at T = 6.0, so that the non-monotonic temperature 
behavior in Fig. ^jp. is not an artifact of our fits, but is a 
genuine effect (see Supplementary Material, Fig. SM-2). 
In addition, we found very similar results for the collec- 
tive relaxation time, Tc{z,T), which further shows that 
this non-monotonic behavior does not sensitively depend 
on the considered observable. Thus, these results give us 
direct evidence that the relaxation processes responsible 
for spatial dynamic correlations have a non-monotonic 
temperature behavior. To our knowledge, all previous 
numerical and experimental studies of spatially heteroge- 
neous dynamics have reported spatial correlations which 
grow as the temperature is decreased and dynamics slows 
down?"-'-^* . 

In Fig. [3] we summarize the temperature dependence 
of the static and dynamic length scales identified above. 
The static length scales show a modest but steady and 
monotonic growth with decreasing temperature, which 
seems to become more pronounced below T ss 6.0. These 
are natural findings from the RFOT theory perspec- 
tive^i^i^, in the sense that static correlations should 
only become prominent below Tc- We also include in 
Fig. [3] the dynamic length scales, which have a strik- 
ing local maximum near T — 6.0. A comparison with 
Fig. [2^ strongly suggests to interpret this maximum in 
the context of RFOT theory in which a dynamic critical 
point at Tc exists in the mean-field limit (but which is 
avoided in finite dimensions) whereas activated dynam- 



ics governed by growing static correlations appears at low 
temperatures. Using RFOT theory Stevenson et al. have 
suggested that around Tc the cooperative domains should 
indeed change shape, in that they have an open structure 
above Tc and a more compact structure below Tc^^ . Our 
findings may be viewed as a striking confirmation of this 
scenario. Although a gradual change from 'flow-like' to 
'hopping' motion was often invoked in the past^, mainly 
to rationalize successes and failures of the mode-coupling 
theory, our results provide a very direct, microscopic evi- 
dence of a change of relaxation mechanism which involves 
no theoretical fitting or indirect interpretation. 

Although the presented results have been obtained for 
a system in the presence of pinned disorder, we empha- 
size that they refiect the behavior of the liquid in the bulk, 
since the nature of the boundary condition does not affect 
averaged static properties. No non-monotonic dynamic 
length scale has been detected so far in bulk systems 
(leaving out systems that show also anomalies in their 
thermodynamics), but we think that this is only related 
to the fact that previous measurements (such as four 
point correlation function x^^"^) only provide a coarse 
representation of spatial correlations, since they mainly 
probe the total number of 'fast relaxing particles' and not 
the details of the shape of the relaxing regions^. There- 
fore the present results should not be taken as a contra- 
diction to previous experiments and simulations, but as 
a new insight into the relaxation mechanism. We hope 
that in the future new experimental techniques will allow 
to detect the non-monotonic T— dependence of £,'^^'^ also 
in real molecular systems. 

If we naively extrapolate our results to lower temper- 
atures, we find a temperature below which static corre- 
lations become larger-ranged than dynamic ones, a sit- 
uation which is physically not very meaningful. Thus 
we are led to speculate that at much lower temperatures, 
dynamic length scales should exhibit an upturn, and per- 
haps become slaved to the static ones, as in the Adam- 
Gibbs picture^ and the scaling regime of RFOT theoryii. 
However, studying numerically this final regime is at 
present too difficult. We suggest that experimental work 
is needed to resolve these issues further. Our study also 
suggests that investigations of confined systems should 
be revisited in both simulations and experiments, and 
the mode-coupling crossover studied more extensively in 
glassformers with different fragility. However, since the 
model investigated here has no unusual features regard- 
ing the relaxation dynamics^lj we expect our results to 
apply also to other simple models such as hard spheres, 
Lennard- Jones-like systems, or soft spheres. 

Methods— We study an equimolar binary mixture of 
harmonic spheres^"* with diameter ratio 1.4, and interac- 
tions between particle i and j given by 



V^J{\n~r,\)^-{l~\n^rJ\yalf if \f,~fj\<<j,, , 

(5) 
where an = ct is the unit of length, cri2 = 1.2, and 



(T22=l-4. The total number of particles is 4320 and all 
of them have the same mass m. Time is expressed in 
units of ^/ma^/e and temperature in units of 10~''£, set- 
ting the Boltzmann constant ks — 1-0. We have used a 
rectangular box of size L^ — Ly — 13.68 and L^ — 34.2, 
yielding a number density p = 0.6749715. This system 
size is sufficiently large to avoid finite size effects. The 
equations of motion have been integrated with the ve- 
locity form of the Verlet algorithm. The longest runs 
extended over 830 million time steps, which took about 
6 weeks of CPU time on a high end processor. In order to 
improve the statistic of the results we have averaged over 



10-30 independent walls. The total amount of computer 
time to obtain the described results was therefore around 
7 years. In practice, to prevent particles to penetrate in 
the frozen half-space, we introduced at the two surfaces of 
the frozen slice an infinitely hard wall, and have checked 
that this frozen geometry has negligible influence on the 
structure of the fluid. (See also Supplementary Mate- 
rial.) Note that we are investigating here a dynamical 
equilibrium, i.e. all fluid particles can leave their initial 
positions and explore the whole confined space, thus as- 
suring thermodynamic equilibrium conditions. 
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Supplementary Material 

Influence of the amorphous waU on the local 
structure of the system 

In order to understand better how the presence of the 
amorphous wall influences the local structure of the con- 
fined liquid we show in Fig.-SM-I a slice of the system 
orthogonal to the confining walls. The two panels corre- 
spond to the two temperature T = 30 and T = 4.5. 

This figure was obtained by making a run that is about 
100 times longer than the a— relaxation time in the bulk 
{=T^^^^{T)) and by superimposing snapshots of the lo- 
cation of the particles every r]*""^. In order to avoid 
overcrowding, only the particles in a slice of thickness 
Ax ~ 1.0 are shown. From this figure one can see that 
in the center of the box the snapshots start to fill up the 
space, i.e. the density of the particles becomes uniform, 
as it should be in a bulk liquid. In contrast to this, the 
presence of the walls influences the local density. How- 
ever, as one can demonstrate analytically, see Ref. [27] 
the structure averaged parallel to the wall is independent 
of z, i.e. is the one of the bulk if one has i) either an 
infinitely extended wall, or ii) makes an average over suf- 
ficiently many independent realizations of the wall. 

Comparing panels a) and b) we see that the range 
over which one finds an influence of the amorphous 
walls on the local structure of the confined liquid does 
indeed grow with decreasing temperature, and in the 
manuscript we demonstrate that this infiuence decreases 
exponentially with z with a temperature dependent 
prefactor (figure lb). To quantify the distance over 
which the wall (locally) modifies the density profile, it 
is reasonable to use ^^tat-mt^ instead of ^s***'^ since the 
former length takes into account also the prefactor of 
the mentioned exponential. 



Time dependence of the overlap at flxed dis- 
tance from the wall 

In order to show that the non-monotonic temperature 
dependence of the dynamic length scale ^s^" is not 
just an effect of the way we have extracted this length 
scale from the data, and that the decrease is present at 
several temperatures, we show in Fig. SM-2 the time 
dependence of the self-overlap qs{z, t) at fixed value of z 
and several temperatures. One recognizes that at high 
temperatures, see curves for T = 10, the curves for the 
different values of z almost coincide, which implies that 
the length scale over which the dynamics is influenced 
by the amorphous walls is very small. With decreasing 
temperature, curves for T — 6.0, one sees a significant 
spread of the curves in the a— regime, indicating that 
^f^" has increased. It the temperature is lowered even 
further (T = 5.25 and T — 5.0) this spread is less 
pronounced than at T = 6.0, i.e. the ^^^^ has decreased 
again. From this figure we thus can conclude that the 
non-monotonic T— dependence of ^f-*'" is not an artifact 
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FIG. 4: Figure SM-1: Superposition of snapsliots of the sys- 
tem at T = 30.0 and T = 4.5 (panel (a) and (b), respectively). 
The vertical dashed lines show the location of the hard walls. 
Only particles in the slab < Xi < 1.0 are shown. Individual 
snapshots are separated by about r]'""^(r) and the total time 
is about 100 r]'"""(r). 



of the way we have analyzed the data and that it is 
seen directly from the time correlation functions at all 
temperatures below T = 6.0. 



Determination of the critical temperature of 
mode-coupling theory 

In Fig. 2a of the paper we show an Arrhenius plot of 
the relaxation times for the system of harmonic spheres. 
Also included in this graph are fits with the predicted 
power-laws of mode-coupling theory [28] . Following stan- 
dard practice in this field, these fits have been obtained in 
the following way: One assumes a value of Tc and makes 
a log-log plot of r]'""^ vs. T — Tc, using of course only 
data point with T > Tc. 
that 



Mode-coupling theory predicts 



bulk 



oc (T - TcY 



(6) 



where the exponent 7 is larger than 1.5, and has for 
hard-sphere-like systems usually a value between 2.0 and 
3.0 [8]. By choosing the value of Tc in such a way as to 
rectify the t^^^^(T) data in a large range, one can de- 
termine Tc- In Fig. SM-3 we show the result of such an 
adjustment for the large particles and we see that the 
power-law describes the data over about three decades. 
We find that for our system the exponent 7 is 2.8, i.e. 




FIG. 5: Figure SM-2: Time dependence of the self-overlap 
qs{z,t) for selected values of z. The temperatures are T — 
10.0, 6.0, 5.25, and 5.0. The curves for T = 5.0 have been 
shifted to the right by a factor of 5.0 in order to avoid over- 
crowding. Note that in the a— relaxation regime the spread 
of the curves is larger for T = 6.0 than for the other temper- 
atures, thus showing the non-monotonic temperature depen- 
dence of the length scale ^^''"(r)- 



a value that is very similar to the ones that are found 
for other glass-forming systems. For the sake of compar- 
ison we have included in the graph also the data for the 
binary Lennard- Jones mixture of Kob and Andersen, a 
system which has been found to be described very well by 
mode-coupling theory [28,29]. The comparison between 
the data for the Lennard-Jones system with the one from 
the harmonic spheres shows that the power-law of mode- 
coupling theory describes the relaxation times well over 
the same number of decades in r]'""^. Therefore we can 
conclude that also for the system of harmonic spheres the 
value of Tc can be estimated with good accuracy. 

We also emphasize that the results we presented here 
that concern the self-quantity {qs{z, t) and r]'""^) are con- 
firmed by the ones for the collective quantities {qc{z,t) 
and T^""^), i.e. there is no relevant difference between 
these two observables as far as relaxation times arc con- 
cerned. 
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